% MITBxByFromBz.m -- Auxiliary MATLAB script for computing Bx and By field maps
%    from Bz field map in the Fourier domain. This script is called by QDMdataprocessing.m
%    
%  References: - R.R. Fu, E.A. Lima, M.W.R. Volk, R. Trubko (2020) "High-sensitivity
%  moment magnetometry with the quantum diamond microscope,"  Geochemistry,
%  Geophysics, Geosystems.  
%              - E. A. Lima and B. P. Weiss (2009) "Obtaining vector magnetic field 
%  maps from single-component measurements of geological samples," Journal
%  of Geophysical Research - Solid Earth, 114, B06102.
%
% -----------------------------------------------------------------------------------
%   Matlab code by Eduardo A. Lima - Copyright (C) 2007-2020 MIT Paleomagnetism Lab
% -----------------------------------------------------------------------------------

function [Bx,By] = MITBxByFromBz(Bz,fs)

%   ---------------------------------------------------------------------------------
%   Bz     -> Z component of the magnetic field measured in a regular planar grid
%   fs     -> Sampling frequency in 1/m
%   ---------------------------------------------------------------------------------

  SHOWGRAPHS=0;      % Set this value to 1 to show frequency response plots, or to 0 otherwise

  [SIZEx, SIZEy] = size(Bz);
  N1=SIZEx;
  N2=SIZEy;

  f1 = [0:N1/2 -(N1/2-1):-1]*fs/N1;     % These freq. coordinates match the fft algorithm
  f2 = [0:N2/2 -(N2/2-1):-1]*fs/N2;

  ff1=(-N1/2:N1/2-1)*fs/N1;             % These freq. coordinates are more suitable to visualization
  ff2=(-N2/2:N2/2-1)*fs/N2;
  
  [F2,F1] = meshgrid(f2+1e-30,f1+1e-30);    % Add an epsilon to frequency variables to avoid going exactly
  [FF2,FF1]=meshgrid(ff2+1e-30,ff1+1e-30);  % over a singularity at the origin of the frequency plane
  kx = 2*pi*F1;
  ky = 2*pi*F2;

  kkx = 2*pi*FF1;
  kky = 2*pi*FF2;
  
  k=sqrt(kx.^2+ky.^2);
  kk=sqrt(kkx.^2+kky.^2);

  etx=-1i*kx./k;                         % Calculate the filter frequency response associated
  evx=-1i*kkx./kk;                       % with the x component

  ety=-1i*ky./k;                         % Calculate the filter frequency response associated
  evy=-1i*kky./kk;                       % with the y component
  
  e = fft2(Bz, N1, N2);

  if SHOWGRAPHS
      figure
      surf(kkx,kky,imag(evx));
      figure
      surf(kkx,kky,imag(evy));
  end
  
  Bx=ifft2(e.*etx,N1,N2,'symmetric');               %Compute Bx map from Bz map
  By=ifft2(e.*ety,N1,N2,'symmetric');               %Compute By map from Bz map
  
